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An analysis method is proposed to study the forward-backward (FB) multiplicity fluctuation in 
high-energy nuclear collisions, built on the earlier work of Bzdak and Teaney. The method allows 
the decomposition of the centrality dependence of average multiplicity from the dynamical event- 
by-event (EbyE) fluctuation of multiplicity in pseudorapidity. Application of the method to AMPT 
and HIJING models shows that the long-range component of the FB correlation is captured by a 
few longitudinal harmonics, with the first component driven by the asymmetry in the number of 
participating nucleons in the two colliding nuclei. The higher-order longitudinal harmonics are found 
to be strongly damped in AMPT compare to HIJING, due to weaker short-range correlations as well 
as the final-state effects present in the AMPT model. Two-particle pseudorapidity correlation reveals 
interesting charge-dependent short-range structures that are absent in HIJING model. The proposed 
method opens an avenue to elucidate the particle production mechanism and early time dynamics in 
heavy-ion collisions. Future analysis directions and prospects of using the pseudorapidity correlation 
function to understand the centrality bias in p+p, p+A and A+A collisions are discussed. 

PACS numbers: 25.75.Dw 


I. INTRODUCTION 

Heavy-ion collisions at RHIC and LHC have two defining characteristics which are the focus of many studies: 1) 
large density fluctuations in the initial state of the collisions that varies event to event, and 2) the rapid formation 
of a strongly coupled quark gluon plasma that expands hydrodynamically with very low specific viscosity. The latter 
characteristic leads to a very efficient transfer of the initial density fluctuations into the final-state collective flow 
correlations in momentum space. Conversely, experimental measurements of the these correlations provide a window 
into the the space-time picture of the collective expansion as well as the medium properties that drives the expansion. 
The measurement of harmonic flow coefficients v n [1-4] and their event-by-event (EbyE) fluctuations [5-7] has placed 
important constraints on the shear viscosity and density fluctuations in the initial state [8-11]. 

Recently, similar ideas have been proposed to study the initial state density fluctuations in the longitudinal direc¬ 
tion [12-15]. These longitudinal fluctuations directly seed the entropy production at very early time of the collisions, 
well before the onset of the collective flow, and appear as correlations of the multiplicity of produced particles sep¬ 
arated in rapidity. For example, EbyE difference between the number of nucleon participants in the target and the 
projectile, iVp art and fVp art may result in a long-range asymmetry of the fireball [13, 14, 16]; the fluctuation of emission 
profile among participants may lead to higher-order shape fluctuations in rapidity [13, 17] (assuming that the emission 
sources for particle production can be associated with individual wounded nucleons). On the other hand, short-range 
correlations can also be generated dynamically including resonance decay, jet fragmentation and Bose-Einstein cor¬ 
relations. These correlations are typically localized over a smaller range of the ry and can be sensitive to final-state 
effects. The longitudinal multiplicity fluctuations, when coupled with the collective transverse expansion, also lead to 
rapidity-dependent EbyE fluctuations of magnitude and the phase of harmonic flow [12, 14, 18, 19]. 

Most previous studies of the longitudinal multiplicity correlation are limited to two rapidity windows symmetric 
around the center-of-mass of the collision system, commonly known as forward-backward (FB) correlations [20, 21]. 
They have been measured experimentally in e + e~ [22], p + p [23-26], p + p [27] and A+A [28, 29] collisions where 
significant FB asymmetric component has been identified. Recently, Refs. [13, 15] generalized the study of the shape 
of the rapidity fluctuation by decompose it into Chebyshev polynomials or into principle components, with each mode 
representing the different components of the measured FB correlation. In this paper, we propose a single-particle 
method that obtains these shape components directly from each event, as well as a two-particle correlation method that 
gives the ensemble RMS-average of these shape components. We apply the method to HIJING [30] and AMPT [31] 
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models and successfully extract the different shape components of the multiplicity fluctuation. The first component 
is found to be directly related to the long-range asymmetry of the fireball, while the higher-order components are 
more related to the short-range correlations. The extracted components are also found to be dampened by the final- 
state interactions. Therefore our method can be used for systematic study of the longitudinal dynamics in heavy-ion 
collisions. 

The structure of the paper is as follows. The next section introduces the method and relates to previous observables. 
Sections III and IV show the properties of the longitudinal shape components extracted from HIJING and AMPT 
models. The meaning of the first few components are discussed within the context of a simple wounded-nucleon and 
particle emission model, and their relations to initial density fluctuations are clarified. Section V compares between 
the single-particle and correlation methods, and a procedure is introduced to further decouple residual centrality 
dependence from the dynamical FB correlations in the correlation function. Section VI discusses new analyses enabled 
by the method, as well as its potential application for understanding the centrality bias effects. 


II. THE METHOD 


The FB correlation can be quantified by two-particle correlation (2PC) function, see for example Ref. [21]: 

(V(? 7 i)V(? 72 )) - (iV(r7i)> 5(771 - 772) 


C{y i,»?2) = 


(N(m))(N( m )) 


(1) 


where the N(rj) = dN/dr/ is multiplicity density distribution in pseudorapidity in one event, the average is over the 
event ensemble, e.g. events within a given centrality class. In experimental analysis, correlation function is usually 
normalized to have an average value of one. The second term in the numerator explicitly removes the self-correlation 
contribution, i.e. one should not correlate a particle with itself. This term is usually dropped in the standard notation, 
since condition 771 + 772 is implicitly assumed, but it is important in our discussion for reasons that will be given below. 

The correlation function can be related to single-particle distribution: 

Citium) = {R(m)R(m ))- , R(v) = 7^4 ( 2 ) 


where R(jf) is the observed multiplicity density distribution in one event normalized by the ensemble average. In the 
absence of EbyE fluctuations, R(r/) = 1 and (7=1. 

One key step in our method is to decompose R(jf) into orthogonal polynomials in the rapidity range [-V,V]: 


<x> / V 

Rid) = 1 + X! a " bS T n (r])i T n (r]) = yjn+ -P n (r)/Y) (3) 

where the Pq(x) = 1 , Pi(x) = x, P2(x) = l/ 2 ( 3 x 2 - 1)..., are Legendre polynomials, and Y characterizes the range 
of the rapidity fluctuations and is chosen to be Y - 6 in current study. The superscript “obs” is used to explicitly 
denote the observed quantity in a single event. The new bases T„(x ) are chosen such that their orthogonality and 
completeness relations are normalized as: 

l / Y / . T n(v) T m(r])dr) = S nm , 1 /Y £ T n (rft)T n (rj 2 ) = 6(771 - V2) ( 4 ) 

n=0 

Our approach is similar to that of Ref. [13] except for two differences: 1) the decomposition is performed on deviation 
from average profile obtained in narrow centrality interval, instead of obtaining (JV(? 7 )} by averging over events with 
different a n values, and 2) the orthogonal bases are Legendre instead of Chebychev polynomials, the latter has a 
weight factor of 1 /\/(1 - (y/Y) 2 ) in the normalization relation that diverges at 77 = ±Y. 

The R(ri) observable provides a natural way to separate the centrality dependence of the (IV( 77 )) from the dynamical 
shape fluctuations for events within fixed centrality: the probability distribution of the N(rj) of all events, p{N(rj)}, 
can be expressed as the sum of the product of the average shape ( N(p)) k and the probability distribution of multiplicity 
shape p{R(r))k} for centrality class “k”: 


p{N(ri)} = Ej; {N(rj)) k p{R(rj)k}■ 


( 5 ) 


Events are first divided into narrow centrality classes according to their total multiplicity M in |? 7 | < Y. Next, the 
average multiplicity distribution (N(rj)) is calculated for each event class, which is then used to calculate the EbyE 
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R(rj). The coefficients of T n and their statistical uncertainty are calculated as: 


obs 


= - (L, 6a. 


obs 



TniVi) 


( 6 ) 


where the sum is over all particles in the event, and 6 n = 1 for n = 0 and 0 otherwise. The <5a° bs characterizes the 
statistical fluctuations due to finite number of particles in the events, and so in principle it can be used to unfold 
the statistical smearing effects in a° bs . In this paper, however, a more robust data-driven method is used to account 
for the smearing of a° bs due to finite number effect: for each real event, a random event is generated with same M 
by sampling the (iV(?;)} and its coefficients a™ are calculated using Eq. 5, which contain only the statistical effects. 
This method provides a simple but self-consistent treatment of the experimental effects. 

Note that, the T n (rf) bases are oscillating functions in pseudorapidity, in a way similar to the azimuthal flow 
harmonics, hence they are referred to as longitudinal harmonics. The non-statistical component of these longitudinal 
harmonics can be obtained after averaging over many events as: 

{a n a m ) = {a n a m ) - ( a n a m ) (7) 

A special case is the diagonal terms: 

( a n) = (« bS ) 2 )-(« n ) 2 )- (8) 

The a n coefficients can also be obtained from the two-particle correlation function: 


2) = 1 + (R«m)R«r ]2 ))-(R ian (yi)R ran (v2)) 

00 

= 1+ E {K s a°^)-(araT))T n (m)T n (m) 

n,m =0 

00 

= 1+ E ( a na m )T rl (77 1 )T m (77 2 ) 

n,m=0 

1 V' / \ T n(Vi)Tm(V 2 ) + T n (ri2)T m (r]i) 

— 1 + ) 

n,m=0 ^ 


( 9 ) 


where we have used the fact that the i? ran (? 7 i) and i? lan (r? 2 ) are uncorrelated except at 771 = 772 - In other words, one 
could construct a correlation function from random events, then it can be shown that: 

c ian (m,V2) ^ (R™(m)R ian (v2)) = 1 + (io) 

This means that the correlation function excluding self-pairs gives directly the ( a n a m ) as the statistical effects drop 
out after averaging pairs over many events. The last part of Eq. 9 is required by < 7 ( 771 , 772 ) = ( 7 ( 772 , 771 ). Furthermore, 
symmetric collision systems such as Pb+Pb require < 7 ( 771 , 772 ) = <7(— 771 , — 772 ), leading to (a n a n+ 1 ) = 0, i. e. odd 
and even harmonics are uncorrelated. The remaining coefficients can be calculated analytically from the correlation 
function as: 

/ x 1 r \nt \ , 1 T n (7?1 ) Trn( 1 T 2 ) + T n ( 7)2 ) T rn (7?-, ) 

(a„a m }=— J [< 7 ( 771 , 772 ) - 1J---^ 771^772 ( 11 ) 

Figure 1 shows the expected shape of the bases in the correlations function, they are plotted assuming ( a n a m ) = 0.01. 
The base for the first term (aiai) is proportional to 771772 and is characterized by quadratic shape along 771 = 772 and 
7 ?i = -772 but with opposite sign (see similar discussion in Ref. [13]). The base for ( 0202 ) is characterized by four sharp 
peaks at the four corners of the correlations function and a broader peak around 771 = 772 ~ 0 . 

The single-particle method denoted by Eqs. 3 and 7 and the correlation method denoted by Eqs. 9 and 11 
are mathematically equivalent. The single-particle method calculates a° bs for each event and hence allows direct 
correlation with its initial geometry in model calculations. Furthermore, it also allow study of possible non-Gaussianity 
in the distribution of a n . On the other hand, the correlation method calculates all ( a n a m ) in a single pass, and 
systematic effects from experiments are easier to control (e.g. via mixed events). 

The discussion above can be generalized into correlations of more than three coefficients, such as (a n a m ai). For the 
single-particle method, it just requires a simple extension of the Eq. 7; while multi-particle correlation functions are 
required for the correlation method, e. g. < 7 ( 771 , 772 , 773 ) 1 . This is an interesting avenue that deserves further studies. 


1 The multi-particle correlation function is closely related to the multi-bin correlator proposed in Ref. [32] 
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FIG. 1: The shape of the first few bases associated with ( a n a m ) in the two-particle correlation function. They are plotted 
assuming { a n a m ) = 0.01. 


To demonstrate the robustness and physics potential of the method, we carried out a detailed simulation study 
using the HIJING [30] and AMPT [31] models. The HIJING model combines the lund-string dynamics for soft 
particle production and hard QCD interaction for high px particle production, which naturally contains many sources 
of long-range and short-range correlations. The AMPT model starts from the particles produced by HIJING, breaks 
them into partons (“string-melting”) and runs them though partonic transport. The partons are then recombined 
to form hadrons at freezeout density, which in turn undergo hadronic transport. The partonic transport processes 
generate significant collective flow and was demonstrated to qualitatively describe the harmonic flow v n in p+A and 
A+A collisions 2 . Therefore measuring the longitudinal harmonics a n in HIJING and AMPT models allows us to 
understand how longitudinal multiplicity fluctuations in the early time are affected by the final-state interactions. 

The HIJING and AMPT data used in this study are generated for Pb+Pb collisions at LHC energy of ^/snn = 2.76 
TeV. All stable particles with px > 0.1 GeV/c in the pseudorapidity range of |r?| < Y = 6 are used. In the default setup, 
events are first sorted into narrow event activity classes based on total multiplicity M, i. e. the M of all events in 
each class is required to differ from the average multiplicity of event class by at most 1%. The N(rj) distribution is 
then obtained for each event and the a° bs coefficients are calculated. At the same time, a random event containing 
M particles is generated according to (N(rj)) and the coefficients a™ are obtained. The same classification is also 
used for 2PC method, however the ( a r a m ) are calculated directly via Eq. 11 without using the random events. This 
event classification procedure in obtaining ( N(rj )) allows a separation of the centrality dependence of the shape of the 
N(rj) distribution (controlled by M) from the shape fluctuations for events with the same M. Hence we can get a 
clearer understanding of the dynamic FB multiplicity fluctuations separated from the overall multiplicity fluctuation. 


2 The model simulation is performed with the string-melting mode with a total partonic cross-section of 1.5 mb and strong coupling 
constant of a s =0.33. This setup has been shown to reproduce the experimental p-\ spectra and v n data at RHIC and the LHC. 
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FIG. 2: The distributions of coefficients for longitudinal Legendre polynomials from real events a „ and random events a a an 
for HIJING (top row) and AMPT (bottom row) events with 6 = 8 fm. The panels in each row correspond to results from n = 1 
to n = 5. 


For comparison, ( N(r ?)) is also obtained using event classes based on either lV part or impact parameter 6 , where much 
stronger EbyE fluctuation is expected for R(rj). 

In the following, we discuss the properties of the a n coefficients based on results obtained from the single-particle 
method. However, most of these results can be also obtained with the 2PC method. 


III. PROPERTIES OF LONGITUDINAL HARMONICS FROM THE SINGLE-PARTICLE METHOD 

Figure 2 shows the EbyE distributions of a° bs for events with fixed impact parameter 6=8 fm, and they are 
compared with distributions obtained from random events a^ an . The differences between the two types distributions 
reflect dynamical fluctuations in a° bs . These differences decrease for larger n, and the rate of decrease is much larger in 
AMPT events than in HIJING events. By n = 5, the distribution for AMPT events is consistent with pure statistical 
fluctuation. From these distributions, the (o^) signals are extracted via Eq. 8 and shown as a function of n in Fig. 3. 
Significant values of a n are seen for all harmonics in HIJING events, while they decrease rapidly and are consistent 
with zero for n > 4 in AMPT events. This difference is mainly due to stronger short-range correlations present in 
HIJING events (see Fig. 13), but could also due to strong viscous damping associated with final-state rescatterings in 
the AMPT model. Figure 4 compares the centrality dependence of the ai, 02 and <23 in HIJING and AMPT models. 
The signal strength increases towards more peripheral collisions and the values from AMPT model are consistently 
smaller than those from HIJING in all centrality ranges. 

In order to find out whether the FB multiplicity fluctuation is related to the difference between IV part and IV® art , 
a° bs is correlated directly with A part , defined as: 


A- 


part 




Xirt 


^art + ^art 


( 12 ) 


The results for 6 = 8 fm from HIJING events are shown in Fig. 5 (results for AMPT events are similar). A strong 
positive correlation between a° hs and A part is observed, suggesting that the FB asymmetry in the multiplicity distri¬ 
bution is indeed driven by the asymmetry in the number of participating nucleons in the two colliding nuclei. A weak 
correlation is also observed between ag bs and A part , suggesting that the FB asymmetry caused by A part contains a 
small non-linear odd component. On the other hand, there is no correlation between <Z 2 bs (rapidity even) and A part 
(rapidity odd) as expected. The width of these distributions are partially due to statistical smearing effects in a° bs , 
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FIG. 3 : The a n vs n from HIJING (left) and AMPT(right) events with 6 = 8 fm. 




FIG. 4 : Centrality dependence of ai (left panel), d2 (middle panel) and 03 (right panel) for HIJING and AMPT events. 



AAA 

npart npart npart 

FIG. 5 : Event-by-event correlation between a)) bs and A par t for n = 1 (left panel), n = 2 (middle panel) and n = 3 (right panel) 
from HIJING events with 6=8 fm. 

which can be removed by a 2D unfolding (leave for a future work). 

Figure 6 (a) compares the centrality dependence of \J(a\) and yj (Ap art ). The similarity in their shapes suggest that 
the asymmetry between fVp art and iV® art is primarily responsible for the FB asymmetry in N{rf). Note that the FB 
asymmetry of R{rj) arising from a\ can be estimated as, Ar(t]) rj \J(a\)Ti(ri) = (a \)|. The results in Fig. 6 (a) 
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FIG. 6 : Comparison between sj(a\) and RMS asymmetry in 7V part , \J{A part ) (left panel), as well as between total multiplicity 
fluctuation in terms of ai\r ch / {N c h) and fluctuation of total lV par t (right panel) in HIJING and AMPT models. 


imply \J(a\) « 0.7a J (dlp al . t }, an< i hence Ar( 6 ) = \J\\J{a\) ~ 0.86^(Ap art ). Therefore, the multiplicity fluctuations in 
the very forward (backward) rapidity (± 6 ) are mostly driven by the fluctuations in N palt (N part ). On the other hand, 
the fluctuation of total multiplicity M is expected to be driven mainly by the fluctuation of N part - N palt + N part . 
Given that ai is driven by N paTt - N part , the fluctuation of M should not be independent from fluctuation of a\. 
Figure 6 (b) compares the relative multiplicity fluctuation, ctm/ {AI), with the fluctuation of number of participants 
^Afpart/ (-/Vpart)- Indeed, the two show very similar centrality dependence after applying a constant scale factor. 

The results shown so far are obtained by calculating (7V(?;)) in narrow bins of M. Figure 7 compare these with 
results obtained in narrow slices of N par t or 6 . This comparison is useful because experiments can only measure 
( N(r])) in finite centrality interval for which the overall multiplicity can still have significant fluctuations. Figure 7 
shows that the values of a± and 03 have very weak dependence on the averaging scheme, while 02 has rather strong 
dependence. The latter suggests that a significant component of the ci 2 obtained for binning in N part or b arises from 
the residual centrality dependence in the shape of ( N(r /)). To see how this residual centrality dependence can arise, 
Fig. 8 compares the ( N(rj )} obtained for events in the upper or lower tails of the total multiplicity distribution for 
all events with 6=8 frn. The ratios on the right panel show that the shape of {N(i 7 )) can still vary significantly 
for events with the same impact parameter but different M , and this variation leads to a significant 02 contribution. 
Nevertheless, after removing this residual centrality dependence by binning events in narrow M ranges, a significant 
a 2 signal still remains. This irreducible 02 could reflect strong event-by-event fluctuations in the amount of nuclear 
stopping or shift of the effective center-of-mass of the collisions [33, 34], Similar results are also seen in HIJING events 
(not shown). 


IV. CORRELATING ai WITH SPECTATOR ASYMMETRY 

If the ai coefficient is correlated with the fluctuations of A r part - N palt , then it should be anti-correlated with the 
asymmetry in the number of spectator nucleons Nj pec - N^ ec since: 

^art - ^art = “(^spec ~ ^spec) ■ (13) 

The number of spectator nucleons can be measured using calorimeters placed very close to the beam-line in the 
forward region. For example, the Zero-degree Calorimeters (ZDC) installed in all RHIC and LHC experiments can 
count the number of spectator neutrons, N neu , in each event with rather good precision. Unfortunately, the measured 
neutrons only constitute a small fraction of all spectator nucleons, and hence the correlation between N paTt - N palt 
and FB neutron asymmetry N aeu - N aeu is expected to be very weak. Nevertheless, studying the correlation between 
ai and Nf pec - N^ ec provides an independent and data-driven way for understanding the origin of the FB multiplicity 
correlations. 
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FIG. 7 : Comparison of the a n obtained from three averaging methods, i.e. binning in total multiplicity, A part or impact 
parameter 6, for ( N(rj )) used in Eq. 2 for n = 1 (left panel), n = 2 (middle panel) and n = 3 (right panel). 




FIG. 8: (a) The average multiplicity distributions for events selected in three multiplicity ranges (see insert) and (b) the ratios 
to the all events. All events are generated for AMPT model with 6 = 8 fm. 


Figure 9(a) shows the ALICE measurement of the correlation of the ZDC energy with ZEM (4.8 < \r]\ < 5.7) energy 
in Pb+Pb collisions at ^/snn = 2.76 TeV [35]. The latter has a very strong correlation with the Silicon Pixel Detector 
(SPD) situated in mid-rapidity (|? 7 | < 1.9) as shown by the insert panel. The ZEM signal can be mapped onto the 
iVp a rt assuming em x i? V part , and the ZDC signal is converted to N neu from the expected energy for each spectator 
nucleon of 1.38 TeV: N neu = Ezdc/ 1-38. From this, the correlation between V part and the average number of neutrons 
(iV n eu) is estimated and shown in Fig. 9(b), where the error bars indicate the approximate standard deviations. This 
correlation is then down-scaled by a factor of two in both axes to give the correlation between 7V part and \V aeu ) or 
between fV® art and \V aeu ). However, the error bar is reduced only by a factor of \/2 assuming the sampling of lV aeu 
is independent of V acu once the values of lVp art and Ad art are fixed in each event (hence N^ pec = 208 - _/V pal . t and 
./V s ®e c = 208-_/V part are also fixed). This new distribution is then used to generate the V aeu and N® eu for each HIJING 
or AMPT event based on its _/V part and N pait values. Finally we calculate the correlation between fV aeu - V aeu and 
a° bs . 

The results of this study for AMPT events is summarized in Fig. 10. A clear anti-correlation is seen in mid-central 
and central collisions. However the correlation is positive in peripheral collisions, which reflects the fact that the value 
of lV neu is positively correlated with V part in the peripheral collisions (see Fig. 9(a)). This correlation is very weak, 
a° hs varies by a few percent in the available range of N^ eu - N® eu , but should be measurable in experiments. 
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FIG. 9: a) The correlation of signals in ZDC and ZEM from ALICE experiment, the insert shows the correlation of signals in 
ZEM and SPD. Then number of neutrons are calculated as N neu = Ezdc/1-38. b) The inferred correlation between N neu and 
fVpart used in this paper. 
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FIG. 10: The estimated correlation between a° bs and N^ eu - N^ eu for peripheral (left panel), mid-central (middle panel) and 
central (right panel) Pb+Pb collisions. 


V. ADDITIONAL INSIGHTS FROM TWO-PARTICLE CORRELATION METHOD 

As discussed in Sec. II, a n coefficients can also be calculated from correlation method via Eq. 11. Figure 11 (a) 
shows the correlation function and (a n a m ) values from AMPT events with b - 8 fm. The shape of the correlation 
function already suggests the dominance of the () term (compare with Fig. 1). The coefficients are compared with 
those obtained from the single-particle method via Eq. 7, identical values are observed. This consistency is expected 
since the two methods are mathematically equivalent. A selected set of coefficients are shown in Fig. 11 (b). No 
correlations are observed between the odd and even coefficients as expected for symmetric collision system, while 
small anti-correlations are observed between odd or even terms, i.e. (a n a n+ 2 ) < 0 and (a n a n+ 4 ) < 0 . 

One important practical advantage of the 2PC method is that it provides a natural way to separate the residual 
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FIG. 11: The correlation function (left) and corresponding spectrum (a„a m ) for n,m < 9 (right panel) for AMPT events 
generated with b = 8 fm, where the {N(r/)) is calculated in narrow multiplicity bins. The spectrum are compared with those 
calculated directly from the single-particle method. 


centrality dependence of average shape of from the dynamical shape fluctuations for events with the same 

centrality. Eq. 9 can be rewritten as : 


1 


1 


C(ih,r] 2 ) = 1+ - {a 0 a 0 )+ — E (a 0 a n ) (T n (r] 2 ) + T n (rji)) + E ( a n a m ) 

^ V 2 n=l n,m =1 


T n (m)Tm{m) + T n (r]2)T m (r) i) 


(14) 


The first term (aoao) reflects the multiplicity fluctuation in the given event class, which drops out from the expression 
if is normalized to have a mean value of one (we shall assume that in the following discussion). The 

second term represents residual centrality dependence in the shape of (N(r])). The last term encodes the dynamical 
shape fluctuations for events with fixed centrality, which can be isolated by dividing the correlation function by its 
projections on the 771 and 772 axes: 


Cn(T/i , 7 / 2 ) = 


C P (r h )C p {rj 2 ) 


n , /C(in,ri2)dri2 f C(r}um)dm 

CM = - 2Y - ,c p (i l2 ) = - 


2 Y 


(15) 

(16) 


The new correlation function ensures that any residual centrality dependence is taken out from the measured coeffi¬ 
cients: 


Cn(? 71,772) = !+ E ( a n«m) 


T n (rn)Tm(m) + T n (r)2)T m (rii) 


(17) 


where the new coefficients are: 


^ {O’n.Um) {doQ'n) (fl'OQ'm) • ( 1 ^) 

They differ from the original coefficients by a small term (aoa n ) (aoa m ), representing the contribution from the residual 
centrality dependence. Alternatively, Cn can also be defined as: 

Cn(?7i,?72) = C(j] 1,772) + 1 - C p (rn)C p (r) 2 ) ( 19 ) 


or: 


CN(vum) = C{i 71 , 772 ) + 2 - C p (r]i) - C p {rj2) ■ 


( 20 ) 
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FIG. 12: The correlation function (left), the product of the projections on two axes (middle) and the redefined correlation 
function via Eq. 15 (right panel) for AMPT events generated with 6 = 8 fm. The (N(r))) is calculated using all events. The 
shallow dip structure shown in the right panel is already present in the left panel. 



FIG. 13: The correlation function defined via Eq. 15 for AMPT (left) and HIJING (right) events generated with 6 = 8 fm. 
The (N(r/)) is calculated using all events. 


Equation 19 practically gives the same answer as Eq. 15. Eq. 20 is not preferred as it does not remove the (aoa n ) (ao<x m ) 
contribution in ( a n a m ), although in practice the relative difference between the two is only a few percent. For all 
results shown below, definition Eq. 15 is used. 

Figure 12 shows the original correlation function, the product of its projections to the two axes, and the renormalized 
correlation function for AMPT events for b = 8 fm, where the average distribution (N(r/)) is calculated in one bin (as 
appose to many narrow multiplicity bins then summed as in Fig. 11). Despite the significant difference in the original 
correlation function due to the residual centrality dependence, the renormalized correlation function is very similar to 
that shown in Fig. 12. The small difference in the four corners of the correlation functions can be attributed to the 
difference in (a^) between different binning schemes shown in Fig. 7(b). Thus the Cn (^ 1 , 772 ) defined in Eq. 17 provides 
a robust way to extract the dynamical shape fluctuations nearly independent of the choice of centrality classes. 

Figure 13 compares the correlation functions between the HIJING and AMPT, the correlation function from AMPT 
appears much broader than the HIJING, which is partially responsible for the faster decrease of the spectrum shown 
in Fig. 3. The AMPT events also show an interesting shallow minimum around Ar/ = 0 with a width of about ±0.4. 
Since it is absent in HIJING events, this structure must reflect the influence of the final-state effects implemented 
in the AMPT model. The correlation function is an intuitive observable for understanding the influence of different 
underlying physics. 
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Note that the correlation function obtained via this procedure is affected by a small bias from short-range component, 
denoted as JsRcCbi, ^ 2 ), via the normalization procedure of Eq. 15. The < 5 src(??i, V 2 ) distribution typically is relatively 
flat along 771+772 with a rather narrow width in the 771-772 direction. In this case, one can easily see that the contribution 
of <5sRc(bij V 2 ) to C p is not uniform in 77 : if the first particle is near mid-rapidity 771 ~ 0 then all pairs in JsRcCffij 772 ) 
contributes to C v (r)\), whereas if the first particle is near the edge of the acceptance 771 ~ ±Y then only half of the 
pairs in SsRciVhVz) contributes to C p (iji). However the short-range component contribution can be estimated, e.g. 
via an experimental procedure discussed in Ref. [36], then such acceptance bias can be removed by redefining the 
projection function and Cn function as: 


^ S ub/ i / [£(771,772) -dsnc(Vi,V 2 )]dri 2 b / [£(771,772) - Ssncirium)] dm 
C p vhJ =-^7-> C P Km) =-^7- 

r „ t x_ C( 771 , 772 ) 

n( )_ C'| Ub (77l)q» Ub (772)‘ 


( 21 ) 

( 22 ) 


Therefore is only corrected for the residual centrality dependence and is free of bias from short-range correlations. 
One can use instead of On to extract a n -spectra. The main effect of the bias is reduce the value of On relative to 
C(j at the four corner’s of the 771,772 phase space. We shall leave this topic for a future study. 


VI. DISCUSSION AND SUMMARY 


We have introduced two complimentary methods for detailed study of the event-by-event fluctuations of particle 
production in the longitudinal direction. The single-particle method gives the coefficients in each event, which can 
be directly relate to the fluctuation of the initial geometry in model calculation. On the other hand, two-particle 
correlation method suppresses the statistical noise on the ensemble basis and hence does not require the construction 
of random events. The correlation method is particularly suitable for small collision system, such as p + p or p+Pb 
collisions, where the EbyE statistical fluctuation is very large. Furthermore, the influence of the detector effects is 
straightforward to remove in the correlation method, and hence it should be considered as the primary method in the 
experimental data analysis. 

The correlation method discussed in this paper can be generalized into correlation between multiplicity of particles 
of any two different types. For example one can measure the correlation between multiplicities for positive and 
negative particles: 


c + (vum) 


(N + ( m )N-(ri 2 )) 

<TV+(7 ?1 ))(iV-(772)>’ 


(23) 


which allow the extraction of (a^a n ). Assuming equal multiplicity for positive and negative particles, the coefficients 
for positive particle a + n and negative particles a“ are related to those for inclusive particles via: 

(«n) = j ((<40 + ( a~a~ ) + 2 {a*a~)) (24) 

Due to local charge conservation effects, the correlation between positive and negative particles is expected to be 
stronger than inclusive correlation. Indeed the AMPT or HIJING simulation studies suggest that {a^a~) > (a^) > 
(<4g 4) = (a“a“). The results shown in Fig. 14 implies that the dip around rji ~ 772 seen in the inclusive correlations for 
AMPT model (e.g. Fig. 11) arises mainly from same-charge pairs, although the opposite-charge pair correlation also 
shows a shallow dip. Such dip is absent in HIJING events independent of the charge combination. These structures 
reflect the important role of the final-state interaction and hardronization mechanism (via simple coalescence in 
AMPT) on the charge-dependent correlations. Note that the charge-dependent correlation function is related to the 
well known balance function B(Arf) [37]: 


2R(At7) = 2C + -(At7) - C ++ (A77) - CT^At?), 


(25) 


The stronger correlation strength for opposite-charge pairs than the same-charge pairs as shown in Fig. 14, implies 
that the balance function should peak around Arj = 771 - 772 = 0 and fall slowly to large Arj (i.e. not sensitive to the 
dips), consistent with earlier observations [38, 39]. 

Similarly, one could also divide particles into high px and low px with equal multiplicity. In this case, the coefficients 
can be written as 

K) ~ \ + («n«n) + 2 (an«n>) 


(26) 
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AMPT b=8fm 
one average bin 

(a) 


opposite charge 



FIG. 14: The correlation functions for same-charge pairs (left panel) and opposite-charge pairs (right panel) for AMPT events 
generated with 6 = 8 fm. 


where and aj) are coefficients for high px and low px particle multiplicity, respectively (for example > 1 GeV/c and < 
1 GeV/c). We observe that (a^a^) > (a^ajj) > (a^a^) (not shown), presumably due to short-range correlations related 
to jet fragmentation, which are stronger for higher px particles. It would be interesting to study the factorization 
behavior of the multiplicity correlation by calculating a factorization ratio, similar to what is often used in azimuthal 
flow correlation analysis [40]: 


a^ah 


\/( a n a n)\/( 


aL a L) 


(27) 


The breaking of the factorization can be used to understand the px dependence of the long-range and short-range 
correlations. 

The a n coefficients can be significantly affected by the short-rangle correlations. One way to suppress such short- 
range correlation is by requiring the pairs to be separated in azimuthal angle </> [21, 26] 3 . However the challenge is 
to understand role of the harmonic flow v n and their EbyE fluctuations, since harmonic flow introduces non-trivial 
multiplicity correlations between particles in different <f> regions. 

In order to study dependence of observables on the size of the collision system, many measurements classify collisions 
according to event activity or centrality in certain r) range. The key challenge in centrality definition is to understand 
dynamical multiplicity correlations between the ?y range used for centrality determination and range used for the 
observable. This is an open issue particularly important in small collision system such as p + p and p+Pb collisions, 
where the bias associated with centrality selection often dominates over the experimental uncertainties [41-45]. Our 
method can be used to measure and quantify such multiplicity correlations, which can then be used to understand the 
influence of centrality biases in other measurements. Since p+Pb is an asymmetric collision system, the correlations 
between odd and even terms may not vanish, which can be studied by measuring (a n a n+ 1 ). 

In summary, a method has been proposed to study the longitudinal multiplicity correlations in high-energy nuclear 
collisions. In this method, events are classified into narrow event activity bins, and EbyE fluctuations are then 
extracted relative to the average multiplicity distribution in each event activity bin. This procedure allows the 
separation of the centrality dependence of the multiplicity distribution from the dynamical shape fluctuations. The 
multiplicity correlations are extracted using the single-particle distribution or two-particle correlation function. The 
extracted signals are decomposed into a set of orthogonal longitudinal harmonics in terms of Legendre polynomials, 
which characterize various components of the multiplicity fluctuation of difference wavelength in i). The first several 


3 In principle, the full information of the transverse and longitudinal multiplicity and flow fluctuations is contained in the 3-D correlation 
function C(?ji,r( 2 , A0). 
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coefficients a n are obtained and found to decrease slowly with n in HIJING model but very rapidly with n in AMPT 
model, which could be due to viscous damping effects of the longitudinal harmonics by the final-state rescattering 
effects. The a\ signal is found to strongly correlated with the asymmetry in the number of forward-going and 
backward-going participating nucleons; while a nonzero <22 signal could be related to the fluctuations of the nuclear 
stopping or shift of the effective center-of-mass of the collisions. This geometrical origin of the ai can be experimentally 
verified by observing an anti-correlation between a\ and the asymmetry of the spectator nucleons detected by the 
zero-degree calorimeters. Two-particle pseudorapidity correlations also reveal interesting charge-dependent short- 
range structures in AMPT model but are absent in HIJING model, suggesting that these structures are sensitive to 
the underlying hadronization mechanism. Hence measurement of the multiplicity fluctuation in terms of longitudinal 
harmonics provide an promising avenue for understanding the particle production mechanism in the early stage of the 
heavy-ion collisions and for probing the final-state rescattering effects. The proposed two-particle correlation method 
is particularly suitable for high-energy proton-lead and proton-proton collisions where the longitudinal multiplicity 
fluctuations are very large and are responsible for the biases in the centrality definition. Since our method correlates 
event activities between separate rapidity ranges, it provides a useful way to unfold and quantify the centrality 
correlations between different rapidity ranges. 

We appreciate fruitful discussions with R. Lacey. This research is supported by NSF under grant number PHY- 
1305037 and by DOE through BNL under contract number DE-SC0012704. 
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